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ABSTRACT 

Iterative methods are considered for the solution of a coupled pair of 
second order elliptic partial differential equations which arise in the field 
of solid state electronics. A finite difference scheme is used which retains 
the conservative form of the differential equations. Numerical solutions are 
obtained in two ways - by multigrid and dynamic alternating direction implicit 
methods. Numerical results are presented which show the multigrid method to 
be an efficient way of solving this problem. 
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Introduction 


In solid state electronics the designers of PIN diodes are Interested in 
the effect on the performance of the diode due to changes in various design 
parameters. An advantage of producing a good mathematical model is that the 
testing, when performed experimentally, could be a lengthy and expensive 
process. A description of the model and the derivation of the equations can 
be found in Altchison and Berz [2] . The model gives rise to a coupled system 
of elliptic partial differential equations. 

In this paper we consider numerical techniques for the solution of this 
pair of equations. We consider a two-dimensional diode which is defined in 
Cartesian co-ordinates and where it is assumed that the diode is very long in 
the third dimension. After deriving the finite difference equations we 
describe applications of a multigrid method and the dynamic A.D.I. (D.A.D.I.) 
method to obtain numerical solutions. 


2. The Differential Equations 

The problem is formulated in terms of the carrier density c(x,y) and a 
stream function u(x,y). The behavior of diodes which are effectively two- 
dimensional and of rectangular cross-section can be described by the following 
equations: 


a 2 a 2 

9 c 3 c 

2 + 2 
9x 9y 


c. 


( 1 ) 


9 (i JLSi') . 1 (I. du ~\ 

3x 'c 3x^ 9y ^c 9y^ 


0 , 


( 2 ) 


in the region R ■ { (x,y) :0<x<A,0<y<B} 
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The boundary conditions are 


3c 

b 

3u 3u 


on 

x = 0, 

(3,4) 

3x ~ 

(1+b) 3y*3x 

3c 

-i 

i 

3 | 

CO \ 

S| 

CO 1 

-(1+b) 3c 

on 

x = A 

(5,6) 

3x 

(1+b) 3y*3x 

b 3y* 

3c _ 
3y 

0, u = 0, 


on 

y = o. 

(7,8) 

3c 

3y 

-sc, u = -1 , 


on 

n 

(9,10) 


where s and b are positive constants. 

There are two quantities in which the designers of diodes are 
particularly interested. The first is the equipotential check, K(y) , given 
by 


K(y) = 


-2b 


( 1+b ) 0 


I ~|^ dx + + l°g(c(0,y) *c(A,y))« (11) 


"b+1 


c(A,y) 


Aitchison [1] showed that this quantity is constant. The second quantity of 
interest is the total charge, Q, defined by 


Q = // cdxdy. 
R 


We can easily verify that Q can be expressed in the form 


- s / c(x,B)dx. 
0 


( 12 ) 


Q - 1 


(13) 
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3* Finite Difference Approximation 

The finite difference equations are constructed using the integration 

method of Varga [8]* This method was also used by Aitchison [1] who solved 

the problem using Newton's method and a sparse matrix routine# The technique 
is used because the conservative form of Eq. (2) is retained in the finite 
difference scheme. 

We consider a rectangle R in which A « 2 and B ■ 4. The region R 

is covered with a square grid of step size h in both the x and y 

directions where Nh « 2. Let c^ j and u^ j be the values of c(x,y) 

and u(x,y) at the grid point where x^ * ih and yj « jh. Let the 

region r^ j be defined as lying within R and being bounded by the lines 


x » x^ - V2 h, x =* : 

Xjl + V 2 h, y = 

yj - V 2 h 

and y = 

yj + V 2 h. 

In 

our 

application of the 

technique the 

regions 

r i,j are 

either square 

or 

rectangular. Various 

regions r i > j 

are shown 

on Fig. 1. 

bet s lfj 

be 

the 


boundary of the region r i j • 

We first consider Eq. (1). Integrating Eq. (1) over the region r^ j 

gives 

IS ( 3 § + ~~~f ~ c)dxdy - 0. (14) 

r i,j 31 3y 

We apply Green's theorem to the first two terms in this integral to obtain 

/ ds - // cdxdy = 0 (15) 

S iJ n r i,j 

where n is the unit outward drawn normal. 

The finite difference approximation at internal points is therefore given 


by 
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- c. 


'1+1 i>J 


)h + (■ 


~ C i,j 


)h+ (■ 


'!~1>J 


- c 


1..1 


)h 



( 16 ) 


where 0<i<N, 0<j < 2N. This equation can be simplified to give 


(17) 


which is the same as that obtained by using the standard five-point finite 
difference approximation to Eq. (1). 



Figure 1 . 
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To construct the finite difference approximation along x ® 0 we need to 
consider the following integral 


| y j+ V 2 l£ . = _b_ r y j+ 1/2 Ju . 

J 9x ay (1+b) J 9y 7 

y j- V 2 V 1/2 


(1+b) 


(u( 0 , (j+ I /2 )h) - u(0, (j- V 2 )h)} 


b f u 0.1+l ~ U 0.1-l j 
(1+b) 1 2 


(IB) 


where 0 < j < 2N. In deriving Eq. (18) we have used boundary condition (3) 
and the approximation u(0 , ( j+ 1/2 )h) = (u Q j+u Q j +1 )/2. Using Eq. (18) we 
obtain the following finite difference approximation to Eq. (1) along x - 0: 

2 c i,j + c o,j+i + c o,j-i ' 4 c o,j " irW K,j+i " u o,j-i^ = h 2 <: o,j • 

where 0 < j < 2N. In a similar fashion we can obtain finite difference 
approximations to Eq . (1) along other parts of the boundary of R. 

We now consider the discretization of Eq. (2). Integrating Eq. (2) over 
the region r^ j yields 

n fj& < 19 > 

r i.J 

Applying Greenes theorem to this equation we obtain 

/ ± & ds - 0. (20) 

9 1,1 

The finite difference approximation to Eq. (2) at internal points is therefore 



6 


U l+_1 J , U i,J+l " U 1,J , U i-l,j ~ U i,l , ~ U l,j 

C i+l,j + C i,j C i,j+1 + C i,j c i-l,J + C i,j C i,j-1 + C i,j 


0, 

( 21 ) 


where 0 < i < N, 0 < j < 2N. 

To construct the finite difference approximation along x - 0 we 
consider the following integral 


/ y J+ 1 /2 

y j-V 2 


1 3u 
c(0,y) 3x 


a“ (0,y)dy =» (1+b)/ 


r J+ V 2 


1 I” (o.y)dy 


c(0,y) 3y 


y j- v 2 

(1+b) {log(c(0, (j+ 1/2 )b) )-log(c (0, ( j- l/ 2 )h) )} 


* (1+b) log 


f. c Q, 1+L_ + 

C 0,j-1 + c 0,j 


(22) 


In this calculation we have made use of boundary condition (4). So the finite 
difference approximation to Eq. (20) along x = 0 is given by 


rff l jl ". “ AIL + J. + (V.1-1 I U Q , 1 ] 1 

L/z K,j + c o,r /2 ^ c o,j+i + c o,j ^ 2 1/2 ^ c o,j-i + C 0,j^ 2 


+ C r 


- (l+b)log(- ^ i- ^ ) = 0, 


0,j-l 


0,j 


where 0 < j < 2N. 

Similarly we obtain finite difference approximations to Eq. (2) along 
x ® 2 using the boundary condition given by Eq. (6). 

Along y = 0 and y ** 4 we have 


u 


1,0 


*= 0, 


u 


i,2N 


- -1, 


and 



7 


respectively where 0 < i < N. 

Let Kj be the discrete form of K((j+V 2 )h) for j - 0,1,..., 2N - 1. 
K((j+V 2 )h) is discretized using the trapezoidal rule. The resulting 

discretization is given by 


-2b 

(1+b) 


2 + 

i-0 12 lc l,j+l C i,j } ° 1 N,j+1 °N,j 


+ lo sK c o,j+l +c 0,j^ C N,j+l +C N,j^’ 


(23) 


where the summation notation is defined by 

A a i + a i + Vi +1/ 2Y 

i«0 

Aitchison [1] shows that Kj is a constant independent of j and so the 
above difference scheme exactly conserves this constant. To calculate the 
total charge Q we discretize Eq. (13), again using the trapezoidal rule. 

A A 

Let Q be the discrete form of Q, then Q is given by 

* N.. 

Q - 1 " sh I c . (24) 

i=0 ’ 


4. A Multigrid Algorithm 

We consider a multigrid method of solution to this coupled system of 
equations using a natural extension of the accommodative Full Approximation 
Storage (F.A.S.) cycling algorithm of Brandt [3]. The multigrid method is a 
numerical strategy to solve partial differential equations by switching 
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between finer and coarser levels of discretization. The characteristic 
feature of the method is the combination of a smoothing step and a coarse grid 
correction. During the smoothing step the residuals are not necessarily 
decreased but smoothed. In the following correction step the discrete 
solution is improved by means of an auxiliary equation on a coarser grid. 
This results in an iterative method that is usually very fast and effective. 
A detailed description of the multigrid method can be found in Brandt [3] and 
Hackbusch [7] . 

Let be a sequence of grids approximating the region 

R = { (x,y) : 0<x<2, 0<y<4} with corresponding mesh sizes h^,...,h M . Let h^ = 

2hj c+ i for k = The problem is discretized on each grid using 

the technique described in the previous section. Let the discrete operators 
k k 

L^ and L^ define the resulting discretizations of equations (1) and (2) 
respectively on G^ where L^ depends on c^. For k < M we solve an 
auxiliary equation on G^ (cf. Algorithm 1). The steps of the algorithm are 
given below. 

Algorithm 1 

(a) Set k = M, the initial working level and choose to be a 

suitable tolerance. Choose initial approximations c^ and u^ to c and 

k k if I, 

u respectively. Set f^ and f^ equal to zero since L^c = 0 and 
L^u^ = 0. On coarser grids fk and f ^ will denote the modified right-hand 
sides (see Eqs . (25) and (26)). 

(b) Set ^ * 10 30 . 

(c) Perform one relaxation sweep over all the equations. Compute the 
J^-norm of the residuals, e^, where 
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, -1 /r. k k k 2 k k k B 2i 

= h /[ HLjC +BL 2 u ” f 2 # 2 


(d) If ^ < e ] c relaxation has sufficiently converged on the 

current level, go to Step (f ) • If not, and If convergence is still fast, i.e. 
e k < w ^ ere D is fixed and chosen later, set e^ » e^ and go to Step 

(c). The parameter n is known as the switching parameter. If convergence 
is slow I.e. e^ > ne^ and we are not on the coarsest grid go to Step (e) • 
If convergence is slow and we are on the coarsest grid go to Step (c) to 
perform another relaxation sweep. 

(e) Decrease k by 1. Transfer the current approximations on level 
k+1 to the new level k as follows : 


c = 


k k+1 

Vi c " 


U = 


_k k+1 

w 


where l£ + ^ denotes some transfer of values from the fine grid. The right- 
hand sides for the new level are defined by 


£ 1 - V + *W f i - L 1 c > 


k+1 , k+1 k+1' 


(25) 


-k T k k /T k f^k+l T k+1 k+1 1 
f 2 -L 2 u +« k+1 (f 2 -t- 2 u ). 


(26) 


The factor 4 appearing in the above equations is a scaling factor which is 

2 

introduced because we multiplied through by h before defining the 
difference operators. Set to be tolerance for the problem on 

the new level where 6 is some parameter. Go to Step (b) . 
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(£) If k “ M, the algorithm is terminated since the problem has been 
solved to the required tolerance. If k < M we correct the approximation on 
the next finer grid . Put 


k+1 

c 

k+1 

u 


k+1 , ..kHz' k _k k+l-\ 

: +I k j* 


k+1 . _k+l/- k _k k+l\ 
U +I k 


where the c^ + ^ ' s and u^ + ^'s on the right-hand sides are the previous 

k+1 

approximations on the level k+1 and 1^ is some interpolation of 
values* Increase k by 1 and go to Step (c). 


Multigrid Components 

We use "non-standard" multigrid techniques introduced by Foerster, Stuben 
and Trottenberg [6] and developed by Foerster and Witsch [5] • 

(i) relaxation 

Pointwise Gauss-Seidel relaxation is used with the points ordered in the 
checkerboard (even-odd) manner* The relaxation of the equations is performed 
in the following order: 

k k k 

(1) relax the equation I^u ■» f ^ at the white (even) points i*e* those 

points (x^y^) for which i + j is even; 

k k k 

(2) relax the equation L^u = f^ at the black (odd) points i*e. those 

points (x^,y ) for which i + j is odd; 

k k k 

(3) relax the equation L^c « f^ at the white points; 

k k k 

(4) relax the equation L^c =» f^ at the black points* 

This is just one of a number of ways of performing checkerboard Gauss-Seidel 
relaxation on these equations. We experimented using several alternatives and 
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found that the above order of relaxation was slightly more efficient than the 

others. At the end of the relaxation sweep the residuals of the equation 
k k k 

LjC ■ fj are zero at the black points since at this stage all the black 
point equations are simultaneously satisfied. 

(ii) fine-to-coarse transfer 

Since checkerboard Gauss-Seidel relaxation produces highly oscillating 

residuals it is not advisable to simply transfer the residuals by injection to 

a coarser grid. Instead we transfer the residuals by full-weighting to the 

coarse grid because the coefficients of Eq. (2) are variable. For this 

k-1 

transfer the operator 1^ is defined by 



k 




l/4v! 


2i,2j 


+1/8 ( V 2i+l,2j +V 2i-l,2j +v 2i,2j-l +V 2i,2j+l) 

+ 1 / 16 ( v 2i+l,2j-l +v 2i-l,2j-l +v 2i+l,2j+l +v 2i-l,2j+l)* 


(iii) coarse- to-f ine transfer 

Bilinear interpolation is used to transfer the correction to the fine 
grid to provide a new approximation there* 


5 • A Dynamic A.D.I* Algorithm 

We now show how the dynamic A*D.I, (D.A.D.I,) method of Doss and Miller 
[4] can be used to obtain a numerical solution to this problem* The A.D*I. 
approach first converts Eqs* (1) and (2) to the parabolic equations 


3c 

3t 


3^c 3^c 

2 + 2 
3x 3y 


c , 


(27) 
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^ (28> 

We assume that as t -*■ 00 the solution of these time-dependent equations tends 
to the steady state solution, if one exists* The parameter A appearing in 
Eq. (28) is used to control the interaction between the equations* When these 
equations are discretized in time it means that, effectively, we use different 
time steps for the two equations. The right-hand sides of Eqs. (27) and (28) 
are discretized using the technique described earlier* We note that the 
systems of equations which we solve in this A*D*I. process are tridiagonal and 
diagonally dominant. This means that the necessary matrix inverses exist and 
the solution of the systems by Gaussian elimination is stable without the need 
for interchanges* Detailed discussions of the A.D.I. method and its 
implementation can be found in Varga [8] and Young [9] • 

Each step of the D. A.D.I. method comprises two double sweeps of the 
A.D.I. iteration with time step At together with a bookkeeping double sweep 
of the A.D.I. iteration with time step 2At. At the end of the step we use a 
computerized strategy to determine how to change At for the next step. Here 
we define a double sweep of the A.D.I. iteration to be a double sweep of 
A.D.I. performed on Eq. (28) followed by a double sweep of A.D.I. performed on 
Eq. (27). The method we have described in this section is a form of decoupled 
D. A.D.I. method where, in a given iteration , we have preferred to iterate 
firstly on u and then on c. A true D. A.D.I. method for this system would 
solve simultaneously for both unknowns. In this case we would solve systems 
which are block tridiagonal with 2x2 blocks. The procedure we use is 
described briefly in Algorithm 2. 
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Algorithm 2 

(a) Choose an Initial time step At = At^. The acceleration parameters 
are then h^/At for Eq. (27) and h^/(XAt) for Eq. (28)* Set k = 0 

where k is the total number of time steps we have advanced. Let e be the 

(k) (k) 

required tolerance. Choose initial approximations c and u • 

(b) Start a step of the D.A.D.I. process with current approximations 

(k) , (k) 

c and u 

(c) Perform two double sweeps of the A.D.I. iteration with time step 

At. Let c (k+4) a nd u (k+4) fo e t ^ e new approximations obtained. Compute 
e, the £ 2 “ norra °f the residuals of the steady state equations. 

(d) If e < £, then the residuals are sufficiently small and the 

* 

algorithm is terminated. If e > e, set At = 2At and determine the 

corresponding acceleration parameters. 

* 

(e) Perform a double sweep of the A.D.I. iteration with time step At 
starting with approximations c^) an d u^) to obtain c (h+4) and u (k+4) 
respectively. This is the bookkeeping part of the D.A.D.I. process. 

(f) Compute the test parameter TP given by 


TP = / [SUM/ASUM] , 


where 


and 


sum = n c (k+4) -c (k+4) a ^ + "u (k+4> -u^ k+4 ^n 2 


ASUM - .c <k+4) -c (k) « 2 2+ «„ <k+4 > -„ (k >. 2 . 


(g) It TP > 0.6, then we reject the present D.A.D.I. step, replace 
At by l/16At and go to Step (b) . If TP < 0.6, then we accept the present 
D.A.D.I. step and change At by the factor of 4,2,/3, V2 > V4 for the next 
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step if TP falls in the intervals (-», 0. 05] , (0.05,0. 1] , (0. 1,0.3] , (0.3, 0.4] , 
(0.4, 0.6] respectively. Go to Step (b) . This is the computerized strategy 
for changing At. 


6. Numerical Results 

When s = 0, Eqs. (1) and (2) together with the boundary conditions (3) 
to (10) possess an analytic solution which is given by 

, s _ (bcosh(2-x)+cosh(x) } 
cU,y) - 4(l+b)sinh(2) ’ 

u(x,y) = - V4Y* 

when A = 2 and B - 4. These functions are used as our initial 

approximation to the solution of the problem for s * 0* Numerical results 
are presented for s = 5. The constant b was given the fixed value 2.7. 

In the multigrid method we define a work unit to be the computational 
work in one relaxation sweep over the finest grid. The step size on the 
coarsest grid is h *= 1. The values of the parameters P and <$ in 

Algorithm 1 were chosen to be 0.5 and 0.3 respectively. It was found that the 
choice of p and <5 was not critical in the sense that values of these 

parameters in the neighborhood of the chosen values produced similar 
efficiency of the algorithm in terms of the number of work units* 

The algorithms were terminated when the °f the residuals was 

less than 10"^. The results in Table I indicate the variation with h of the 

values of c and u at the center of the diode, the values of c at the 

corners, and the values of the constants Q and K. In Table II we give 
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details of the multigrid method of solution. Details of the D.A.D.I. method 
of solution are given in Table III for X = 1 and X = 0.05. The results 

were obtained on the Oxford University ICL 2980 computer. 

Various values of X were tried in the D.A.D.I. method and it was found, 

by experiment, that the value X = 0.05 produced the fastest convergence. It 

can be seen from the computational details in Table III that this value of X 
is considerably better than X = 1 for the smallest mesh size h = 0.125. 

However, even with this value of A, the multigrid method performs much better 
than the D.A.D.I. method on this problem. 

The values of the asymptotic convergence factor are a little higher than 
typical for multigrid methods. Experiments were performed fixing each of the 
variables In turn to determine how the convergence rate behaves for a single 
equation. The results of this investigation are given in Table IV. In this 
table we give the asymptotic convergence factors per work unit for different 
values of h. We see that for the single equations the typical multigrid 
rates are realized. Closer examination revealed that the higher convergence 
factors of the system were due to the coupling through the boundary 
conditions. A more effective treatment of the boundary conditions in the 
multigrid context will be the subject of future work. 

The methods described here are not restricted to use on square grids. 
These grids were chosen since there were no advantages in using non-uniform 
grids for this problem. A qualitative discussion of what one learns about the 
diode from this solution appears in Aitchison [1]. 
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TABLE I 



Details of the solution 

for different values of h. 



h = 0.5 

h - 0.25 

h - 0.125 

/-s 

CM 

a 

0.0890 

0.0906 

0.0910 

CM 

t—i 

V 

-0.5542 

-0.5550 

-0.5520 

c(0,0) 

0.2018 

0.2071 

0.2086 

c(2,0) 

0.1203 

0.1228 

0.1253 

c(0,4) 

0.0591 

0.0682 

0.0749 

c(2,4) 

0.0196 

0.0203 

0.0209 

Q 

0.7831 

0.7871 

0.7884 

K 

-1.6415 

-1.6074 

-1.5984 


TABLE II 


Details of Multigrid Method 

h 

Number of work units 

Time (secs.) 

0.5 

46 

0.9 

0.25 

56 

2.3 


0.125 


65 


7.3 
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TABLE III 


Details of D.A.D.I. Method 

h 

Number 

of D.A.D.I. steps 

Time 

(secs .) 



X = 1 

X = 0.05 

X = 1 

X = 0.05 


0.5 

84 

64 

2.3 

1.9 


0.25 

140 

98 

14.0 

9.7 


0.125 

216 

122 

155.7 

88.6 



TABLE IV 



Asymptotic 

convergence factors 



h = 0.5 

h - 0.25 

h = 0.125 

c equation 

0.42 

0.45 

0.46 

u equation 

0.45 

0.35 

0.35 

system 

0.72 

0.75 

0.77 
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